https://github.com/cran/VineCopula - Copula family
setwd('~/Documents/VaR')
prices <- read.csv('./data/StockIndexData.csv')
library(VineCopula)
library(copula)
##
## Attaching package: 'copula'
## The following object is masked from 'package:VineCopula':
##
## fitCopula
library(tseries)
prices.DJIA <- prices[,c('DJIA')]
summary(prices.DJIA)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6547 10220 11020 12020 13270 18640
n <- length(prices.DJIA)
rtn.DJIA <- rep(0,n-1)
for(i in seq(n-1)){
rtn.DJIA[i] <- log(prices.DJIA[i+1]/prices.DJIA[i]) * 100.
}
prices.GSPC <- prices[,c('GSPC')]
summary(prices.GSPC)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 676.5 1133.0 1296.0 1364.0 1478.0 2190.0
n <- length(prices.GSPC)
rtn.GSPC <- rep(0,n-1)
for(i in seq(n-1)){
rtn.GSPC[i] <- log(prices.GSPC[i+1]/prices.GSPC[i]) * 100.
}
prices.NDX <- prices[,c('NDX')]
summary(prices.NDX)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 938.5 1569.0 1944.0 2394.0 2954.0 4910.0 1063
n <- length(prices.NDX)
rtn.NDX <- rep(0,n-1)
for(i in seq(n-1)){
rtn.NDX[i] <- log(prices.NDX[i+1]/prices.NDX[i]) * 100.
}
prices.GDAXI <- prices[,c('GDAXI')]
summary(prices.GDAXI)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2203 4934 6172 6461 7583 12370 60
n <- length(prices.GDAXI)
rtn.GDAXI <- rep(0,n-1)
for(i in seq(n-1)){
rtn.GDAXI[i] <- log(prices.GDAXI[i+1]/prices.GDAXI[i]) * 100.
}
prices.FCHI <- prices[,c('FCHI')]
summary(prices.FCHI)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2403 3649 4220 4291 4836 6857 50
n <- length(prices.FCHI)
rtn.FCHI <- rep(0,n-1)
for(i in seq(n-1)){
rtn.FCHI[i] <- log(prices.FCHI[i+1]/prices.FCHI[i]) * 100.
}
prices.SSEC <- prices[,c('SSEC')]
summary(prices.SSEC)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1011 1572 2104 2288 2850 6092 317
n <- length(prices.SSEC)
rtn.SSEC <- rep(0,n-1)
for(i in seq(n-1)){
rtn.SSEC[i] <- log(prices.SSEC[i+1]/prices.SSEC[i]) * 100.
}
prices.SENSEX <- prices[,c('SENSEX')]
summary(prices.SENSEX)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2600 4757 13670 13040 18840 29680 202
n <- length(prices.SENSEX)
rtn.SENSEX <- rep(0,n-1)
for(i in seq(n-1)){
rtn.SENSEX[i] <- log(prices.SENSEX[i+1]/prices.SENSEX[i]) * 100.
}
rtn.DJIA.ar <- ar(rtn.DJIA, order.max = 1, method = 'ols')
rtn.DJIA.ar.resid.garch <- garch(rtn.DJIA.ar$resid[2:length(rtn.DJIA.ar$resid)],trace=FALSE)
rtn.GSPC.ar <- ar(rtn.GSPC, order.max = 1, method = 'ols')
rtn.GSPC.ar.resid.garch <- garch(rtn.GSPC.ar$resid[2:length(rtn.GSPC.ar$resid)],trace=FALSE)
rtn.NDX.1<-na.omit(rtn.NDX)
rtn.NDX.ar <- ar(rtn.NDX.1, order.max = 1, method = 'ols')
rtn.NDX.ar.resid.garch <- garch(rtn.NDX.ar$resid[2:length(rtn.NDX.ar$resid)],trace=FALSE)
rtn.GDAXI.1<-na.omit(rtn.GDAXI)
rtn.GDAXI.ar <- ar(rtn.GDAXI.1, order.max = 1, method = 'ols')
rtn.GDAXI.ar.resid.garch <- garch(rtn.GDAXI.ar$resid[2:length(rtn.GDAXI.ar$resid)],trace=FALSE)
rtn.FCHI.1<-na.omit(rtn.FCHI)
rtn.FCHI.ar <- ar(rtn.FCHI.1, order.max = 1, method = 'ols')
rtn.FCHI.ar.resid.garch <- garch(rtn.FCHI.ar$resid[2:length(rtn.FCHI.ar$resid)],trace=FALSE)
rtn.SSEC.1<-na.omit(rtn.SSEC)
rtn.SSEC.ar <- ar(rtn.SSEC.1, order.max = 1, method = 'ols')
rtn.SSEC.ar.resid.garch <- garch(rtn.SSEC.ar$resid[2:length(rtn.SSEC.ar$resid)],trace=FALSE)
rtn.SENSEX.1<-na.omit(rtn.SENSEX)
rtn.SENSEX.ar <- ar(rtn.SENSEX.1, order.max = 1, method = 'ols')
rtn.SENSEX.ar.resid.garch <- garch(rtn.SENSEX.ar$resid[2:length(rtn.SENSEX.ar$resid)],trace=FALSE)
rm(prices.DJIA, prices.GSPC, prices.NDX, prices.GDAXI, prices.FCHI,prices.SSEC, prices.SENSEX)
rm(rtn.DJIA, rtn.GSPC, rtn.NDX, rtn.GDAXI, rtn.FCHI,rtn.SSEC, rtn.SENSEX)
rm(rtn.NDX.1, rtn.GDAXI.1, rtn.FCHI.1,rtn.SSEC.1, rtn.SENSEX.1)
Standardize the residual from AR(1)XGARCH(1,1)
std.residual.DJIA <- rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.DJIA.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GSPC <- rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GSPC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.NDX <- rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.NDX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GDAXI <- rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GDAXI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.FCHI <- rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.FCHI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SSEC <- rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SSEC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SENSEX <- rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SENSEX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch
## $fitted.values): longer object length is not a multiple of shorter object
## length
rtns <- cbind(std.residual.DJIA[1067:4532],
std.residual.GSPC[1067:4532],
std.residual.NDX[2:3467],
std.residual.GDAXI[953:4418],
std.residual.FCHI[967:4432],
std.residual.SSEC[658:4123],
std.residual.SENSEX[674:4139]
)
rm(std.residual.DJIA,std.residual.GSPC,std.residual.NDX,std.residual.GDAXI ,std.residual.FCHI,
std.residual.SSEC,std.residual.SENSEX)
cor(rtns,method='pearson')
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00000000 0.97482107 0.4626361048 -0.02897670 -0.0448450319
## [2,] 0.97482107 1.00000000 0.4877918823 -0.02413000 -0.0454992144
## [3,] 0.46263610 0.48779188 1.0000000000 -0.05170464 0.0005284583
## [4,] -0.02897670 -0.02413000 -0.0517046389 1.00000000 -0.0142933609
## [5,] -0.04484503 -0.04549921 0.0005284583 -0.01429336 1.0000000000
## [6,] -0.02312062 -0.02512210 0.0025580560 0.01929367 0.0158593348
## [7,] -0.01925440 -0.02025158 -0.0048500626 0.01343810 -0.0123316220
## [,6] [,7]
## [1,] -0.023120619 -0.019254398
## [2,] -0.025122102 -0.020251582
## [3,] 0.002558056 -0.004850063
## [4,] 0.019293672 0.013438099
## [5,] 0.015859335 -0.012331622
## [6,] 1.000000000 0.013798025
## [7,] 0.013798025 1.000000000
cor(rtns,method='kendall')
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000000 0.840936111 0.339267541 -0.018745030 -0.020497615
## [2,] 0.840936111 1.000000000 0.370060343 -0.012876103 -0.020819688
## [3,] 0.339267541 0.370060343 1.000000000 -0.026871801 -0.001039327
## [4,] -0.018745030 -0.012876103 -0.026871801 1.000000000 -0.012353858
## [5,] -0.020497615 -0.020819688 -0.001039327 -0.012353858 1.000000000
## [6,] -0.003618578 -0.004117842 -0.000861471 0.010674214 -0.002154094
## [7,] -0.004118841 -0.003936654 -0.004190783 -0.001609867 -0.020418012
## [,6] [,7]
## [1,] -0.003618578 -0.004118841
## [2,] -0.004117842 -0.003936654
## [3,] -0.000861471 -0.004190783
## [4,] 0.010674214 -0.001609867
## [5,] -0.002154094 -0.020418012
## [6,] 1.000000000 0.009381258
## [7,] 0.009381258 1.000000000
cor(rtns,method='spearman')
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000000 0.962107436 0.454675399 -0.027960277 -0.030258828
## [2,] 0.962107436 1.000000000 0.484046514 -0.019486288 -0.030636594
## [3,] 0.454675399 0.484046514 1.000000000 -0.040064311 -0.001568855
## [4,] -0.027960277 -0.019486288 -0.040064311 1.000000000 -0.018482587
## [5,] -0.030258828 -0.030636594 -0.001568855 -0.018482587 1.000000000
## [6,] -0.005895352 -0.006491070 -0.001236738 0.015962569 -0.002794678
## [7,] -0.006206004 -0.006164209 -0.006189422 -0.002249416 -0.030567593
## [,6] [,7]
## [1,] -0.005895352 -0.006206004
## [2,] -0.006491070 -0.006164209
## [3,] -0.001236738 -0.006189422
## [4,] 0.015962569 -0.002249416
## [5,] -0.002794678 -0.030567593
## [6,] 1.000000000 0.013779081
## [7,] 0.013779081 1.000000000
splom2(rtns)
g <- rtns[,1]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 1 10 70 366 1183 1401 377 53 3 2
##
## $density
## [1] 0.0002885170 0.0028851702 0.0201961916 0.1055972302 0.3413156376
## [6] 0.4042123485 0.1087709175 0.0152914022 0.0008655511 0.0005770340
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,2]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 4 14 82 362 1139 1409 398 52 4 2
##
## $density
## [1] 0.001154068 0.004039238 0.023658396 0.104443162 0.328620889
## [6] 0.406520485 0.114829775 0.015002885 0.001154068 0.000577034
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,3]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
##
## $counts
## [1] 1 3 13 132 423 1072 1259 464 90 7 1 1
##
## $density
## [1] 0.0002885170 0.0008655511 0.0037507213 0.0380842470 0.1220427005
## [6] 0.3092902481 0.3632429313 0.1338718984 0.0259665320 0.0020196192
## [11] 0.0002885170 0.0002885170
##
## $mids
## [1] -5.5 -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 5.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,4]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 3 19 119 431 1051 1258 478 88 16 3
##
## $density
## [1] 0.0008655511 0.0054818234 0.0343335257 0.1243508367 0.3032313907
## [6] 0.3629544143 0.1379111368 0.0253894980 0.0046162724 0.0008655511
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,5]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 2 20 122 422 1118 1205 479 82 10 6
##
## $density
## [1] 0.000577034 0.005770340 0.035199077 0.121754183 0.322562031
## [6] 0.347663012 0.138199654 0.023658396 0.002885170 0.001731102
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,6]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
##
## $counts
## [1] 1 1 10 42 112 424 1099 1166 460 123 22 4 2
##
## $density
## [1] 0.000288517 0.000288517 0.002885170 0.012117715 0.032313907
## [6] 0.122331218 0.317080208 0.336410848 0.132717830 0.035487594
## [11] 0.006347374 0.001154068 0.000577034
##
## $mids
## [1] -6.5 -5.5 -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 5.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,7]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 3 3 65 373 1228 1346 396 48 3 1
##
## $density
## [1] 0.0008655511 0.0008655511 0.0187536065 0.1076168494 0.3542989036
## [6] 0.3883439123 0.1142527409 0.0138488171 0.0008655511 0.0002885170
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
m <- pobs(as.matrix(cbind(rtns)))
splom2(m)
n.cop <- normalCopula(dim=7, dispstr = 'un')
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
## Estimate Std. Error z value Pr(>|z|)
## rho.1 9.714e-01 6.833e-04 1421.637 < 2e-16 ***
## rho.2 4.661e-01 1.195e-02 38.994 < 2e-16 ***
## rho.3 -3.060e-02 1.702e-02 -1.798 0.07225 .
## rho.4 -3.929e-02 1.701e-02 -2.310 0.02092 *
## rho.5 -1.696e-02 1.705e-02 -0.995 0.31989
## rho.6 -1.430e-02 1.705e-02 -0.839 0.40166
## rho.7 4.915e-01 1.157e-02 42.464 < 2e-16 ***
## rho.8 -2.463e-02 1.703e-02 -1.447 0.14801
## rho.9 -4.021e-02 1.701e-02 -2.364 0.01809 *
## rho.10 -1.877e-02 1.705e-02 -1.101 0.27075
## rho.11 -1.387e-02 1.705e-02 -0.814 0.41593
## rho.12 -4.980e-02 1.699e-02 -2.931 0.00338 **
## rho.13 1.646e-06 1.705e-02 0.000 0.99992
## rho.14 1.580e-03 1.705e-02 0.093 0.92620
## rho.15 -5.006e-03 1.706e-02 -0.293 0.76914
## rho.16 -1.738e-02 1.705e-02 -1.019 0.30805
## rho.17 1.913e-02 1.705e-02 1.122 0.26171
## rho.18 8.482e-03 1.705e-02 0.497 0.61893
## rho.19 9.321e-03 1.705e-02 0.547 0.58467
## rho.20 -1.856e-02 1.705e-02 -1.089 0.27628
## rho.21 1.395e-02 1.705e-02 0.818 0.41315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is 5474
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 196 28
rho <-coef(fit)
copula_dist <- mvdc(copula=normalCopula(rho,dim=7,dispstr = 'un'),
margins=c("norm","norm","norm","norm","norm","norm","norm"),
paramMargins=list(list(mean=mean(rtns[,1]), sd=sd(rtns[,1])),
list(mean=mean(rtns[,2]), sd=sd(rtns[,2])),
list(mean=mean(rtns[,3]), sd=sd(rtns[,3])),
list(mean=mean(rtns[,4]), sd=sd(rtns[,4])),
list(mean=mean(rtns[,5]), sd=sd(rtns[,5])),
list(mean=mean(rtns[,6]), sd=sd(rtns[,6])),
list(mean=mean(rtns[,7]), sd=sd(rtns[,7]))
)
)
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='normal-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,3],main='normal-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,4],main='normal-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,5],main='normal-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,6],main='normal-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,7],main='normal-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
t.cop <- tCopula(dim=7, dispstr = 'un')
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
## Estimate Std. Error z value Pr(>|z|)
## rho.1 0.9713667 0.0007552 1286.156 < 2e-16 ***
## rho.2 0.5145676 0.0124411 41.360 < 2e-16 ***
## rho.3 -0.0293183 0.0177551 -1.651 0.09869 .
## rho.4 -0.0316248 0.0178154 -1.775 0.07588 .
## rho.5 -0.0133095 0.0177248 -0.751 0.45272
## rho.6 -0.0034288 0.0176745 -0.194 0.84618
## rho.7 0.5414400 0.0120706 44.856 < 2e-16 ***
## rho.8 -0.0242672 0.0177692 -1.366 0.17204
## rho.9 -0.0333458 0.0178330 -1.870 0.06150 .
## rho.10 -0.0159670 0.0177309 -0.901 0.36785
## rho.11 -0.0038308 0.0176985 -0.216 0.82864
## rho.12 -0.0473888 0.0177987 -2.662 0.00776 **
## rho.13 -0.0045290 0.0179107 -0.253 0.80037
## rho.14 -0.0062200 0.0177792 -0.350 0.72645
## rho.15 0.0038222 0.0177604 0.215 0.82960
## rho.16 -0.0198613 0.0175692 -1.130 0.25828
## rho.17 0.0191351 0.0175063 1.093 0.27438
## rho.18 0.0013937 0.0175002 0.080 0.93653
## rho.19 0.0018595 0.0175989 0.106 0.91585
## rho.20 -0.0199084 0.0174932 -1.138 0.25509
## rho.21 0.0196046 0.0174675 1.122 0.26172
## df 10.5921609 0.6350647 16.679 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is 5681
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 184 33
rho <- coef(fit)[1:21]
df <- coef(fit)[22]
copula_dist <- mvdc(copula=normalCopula(rho,dim=7,dispstr = 'un'),
margins=c("t","t","t","t","t","t","t"),
paramMargins=list(list(df = df),
list(df = df),
list(df = df),
list(df = df),
list(df = df),
list(df = df),
list(df = df)
)
)
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='t-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,3],main='t-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,4],main='t-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,5],main='t-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,6],main='t-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,7],main='t-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
c.cop <- claytonCopula(dim=7)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
## Estimate Std. Error z value Pr(>|z|)
## param 0.095633 0.005807 16.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is 226.4
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 24 6
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=7),
margins=c("norm","norm","norm","norm","norm","norm","norm"),
paramMargins=list(list(mean=mean(rtns[,1]), sd=sd(rtns[,1])),
list(mean=mean(rtns[,2]), sd=sd(rtns[,2])),
list(mean=mean(rtns[,3]), sd=sd(rtns[,3])),
list(mean=mean(rtns[,4]), sd=sd(rtns[,4])),
list(mean=mean(rtns[,5]), sd=sd(rtns[,5])),
list(mean=mean(rtns[,6]), sd=sd(rtns[,6])),
list(mean=mean(rtns[,7]), sd=sd(rtns[,7]))
)
)
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='Clayton-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,3],main='Clayton-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,4],main='Clayton-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,5],main='Clayton-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,6],main='Clayton-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,7],main='Clayton-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
c.cop <- gumbelCopula(dim=7)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
## Warning in .local(copula, tau, ...): tau is out of the range [0, 1]
fit
## fitCopula() estimation based on 'maximum likelihood'
## and a sample of size 3466.
## Estimate Std. Error z value Pr(>|z|)
## param 1.043675 0.004049 257.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## The maximized loglikelihood is 89.32
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 39 10
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=7),
margins=c("norm","norm","norm","norm","norm","norm","norm"),
paramMargins=list(list(mean=mean(rtns[,1]), sd=sd(rtns[,1])),
list(mean=mean(rtns[,2]), sd=sd(rtns[,2])),
list(mean=mean(rtns[,3]), sd=sd(rtns[,3])),
list(mean=mean(rtns[,4]), sd=sd(rtns[,4])),
list(mean=mean(rtns[,5]), sd=sd(rtns[,5])),
list(mean=mean(rtns[,6]), sd=sd(rtns[,6])),
list(mean=mean(rtns[,7]), sd=sd(rtns[,7]))
)
)
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(rtns[,1],rtns[,2],main='Gumbel-copula DJIA-GSPC')
points(sim[,1],sim[,2],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,3],main='Gumbel-copula DJIA-NDX')
points(sim[,1],sim[,3],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,4],main='Gumbel-copula DJIA-GDAXI')
points(sim[,1],sim[,4],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,5],main='Gumbel-copula DJIA-FCHI')
points(sim[,1],sim[,5],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,6],main='Gumbel-copula DJIA-SSEC')
points(sim[,1],sim[,6],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)
plot(rtns[,1],rtns[,7],main='Gumbel-copula DJIA-SENSEX')
points(sim[,1],sim[,7],col='purple')
legend('bottomright',c('Observed','Simulated'),col=c('black','purple'),pch=21)